(i2> United States Patent 

Chen et al. 



i nil uim i uii in mi idi mil mn nil no inn iiii i m 

US006687422B1 

(io) Patent No,: US 6,687,422 Bl 
(45) Date of Patent: Feb. 3, 2004 



(54) DISCRETE IMAGE DATA INTERPOLATION 
USING SUB-UNITY SHIFT 

(75) Inventors: Qin-Sheng Chen, Beachwood, OH 

(US); Martin S. Weinhous, Moreland 
Hills, OH (US) 

(73) Assignee: The Cleveland Clinic Foundation, 

Cleveland, OH (US) 

( * ) Notice: Subject to any disclaimer, the term of this 
. patent is extended or adjusted under 35 
U.S.C. 154(b) by 0 days. 

(21) Appl.No.: 09/623,781 

(22) PCT Fiied: Feb. 16, 2000 

(86) PCT No.: PCT/USOO/03995 
§371 (c)(1), 

(2), (4) Date: Dec. 1, 2000 

(87) PCT Pub. No.: WO00/49596 
PCT Pub. Date: Aug. 24, 2000 

(51) Int. CI. 7 G06K9/32 

(52) U.S. CI 382/300; 382/280; 345/606; 

708/403 

(58) Field of Search 382/276, 277, 

382/278, 279, 280, 299, 300; 708/403, 
404, 405; 345/606, 615 

(56) References Cited 

U.S. PATENT DOCUMENTS 



5,170,170 A * 12/1992 Soumekh 342/179 

5,481,275 A * 1/1996 Mical et al 345/698 

5,566,284 A * 10/1996 Wakayama 345/587 

5,917,940 A * 6/1999 Okajima et al 382/173 

6,178,271 Bl * 1/2001 Maas, III 382/294 

6,181,135 Bl * 1/2001 Shen 324/309 

6,456,745 Bl * 9/2002 Bruton et al *. 382/298 

6,507,859 Bl * 1/2003 Omori et al 708/300 



OTHER PUBLICATIONS 

"Image Approximation by Variable Knot Bicubic Splines." 
Dennis G. McCaughey, and Harry C. Andrews, IEEE Trans- 
actions on Pattern Analysis and Machine Intelligence, vol. 
PAMI-3, No. 3, May 1981, pp. 299-310. 
"Simple Algorithms and Architectures for B^Spline Inter- 
polation." P.V. Sankar and L.A. Ferrari, IEEE Transactions 
on Pattern Analysis and Machine Intelligence, vol. 10, No. 
2, Mar. 1988, pp. 271-276. 

"Fast B-Spline Transforms for Continuous Image Repre- 
sentation and Interpolation." Michael Unser, Akram 
Aldroubi, Murray Eden, IEEE Transactions on Pattern 
Analysis and Machine Intelligence, vol. 13, No. 3, Mar. 
1991, pp. 277-285. 

"Quadratic Interpolation for Image Resampling." Neil A. 
Dodgson, IEEE Transactions on Image Processing, vol. 6, 
No. 9, Sep., 1997, pp. 1322-1326. 

(List continued on next page.) 

Primary Examiner — Samir Ahmed 

Assistant Examiner — Anand Bhatnagar 

(74) Attorney, Agent, or Firm— Ned Pejic; Calfee, Halter & 

Griswold, LLP 

(57) ABSTRACT 

A method for interpolating digital images is provided. An 
original digital image in the spatial domain is transformed 
into a frequency domain representation via a Fourier trans- 
form. The Fourier transform or frequency domain represen- 
tation is then modified by phase shift terms corresponding to 
image shifts in the spatial domain with sub -unity distances 
matching the locations where the image values need to be 
restored. The original and the shifted images are then 
interspersed together^ yielding an interpolated image. The 
method is particularly useful for two and three-dimensional 
computer tomography (CT) and magnetic resonance imag- 
ing (MRI) images, as well as other medical and non-medical 
digital images. 
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DISCRETE IMAGE DATA INTERPOLATION 
USING SUB-UNITY SHIFT 

FIELD OF THE INVENTION 

The invention relates generally to digital images and 
image processing and, more particularly, to a method of 
. interpolating missing elements from originally sampled 
images and generating an interpolated image therefrom. 

BACKGROUND OF THE INVENTION 

With continuous improvements of imaging and comput- 
ing technologies, more and more medical images are pro- 
vided in digital form. Digital images can be easily manipu- 
lated on a computer to enhance or extract some properties 
which are clinically valuable, but may not be available with 
conventional images on films. In advanced radiotherapy 
treatment planning for instance, the computer tomography 
(CT) number in a voxel is used to determine the tissue 

, attenuation which affects the radiation deposition in the 
tissue. In many clinical applications, such as visualization, 
quantitative measurements, and analysis, it is frequently 
necessary to interpolate an image so that different image 
properties can be explored. 

Image interpolation is a process used to- estimate the 
missing elements from the original sampling grid or image. 
If the sampling rate of a data set satisfies the Nyquist 
sampling criterion, the missing data in principle can be fully 
recovered by convo luting the discrete data with the continu- 
ous interpolating function (sine function). Since the sine 
function is an infinite function while digital images are 

' always confinecf within a finite field, the interpolation can 
only be approximated. A major drawback in image imple- 
mentation using the sine function however is the computa- 
tion efficiency. 

In order to improve the computation efficiency with 
acceptable interpolation quality, different "interpolation 
approaches have been developed. The most popular scheme 
is based on convolution techniques with selected kernels. 
The simplest kernel is the nearest-neighbor function, where 
the values of the interpolating points are copied from the 
nearest samples. Another widely used kernel is the first order 

. polynomial function, where the value at an interpolating 
point is a linear combination of the adjacent samples. A 
major problem with these kernels is the unsatisfactory 
interpolation quality with staircase (i.e., blocky) artifacts in 
the resulting images. More complex kernels in this category 
are the quadratic and cubic spline functions. It has been 
showed that a spline function approximates a truncated 
version of the sine function. Since a spline function is a 
piecewise continuous polynomial function, complexity and 
computing efficiency are the drawbacks. Improper selection 
of the spline function parameters may also significantly blur 
the resulting image. Alternatively, wavelet-based and Gaus- 

* sian kernels have also been proposed. These can be sum- 
marized as attempts to improve the approximation of the 
truncated sine function. Another class of image interpolation 
techniques is based on object connectivity and relative 
displacement. All of the aforementioned processing tech- 
niques suffer from the undesirable effects of long computa- 
tion times and computational inefficiency that make these 
approaches impractical in routine clinical applications. 
Consequently, an image process method that does not suffer 
from such disadvantages is highly desirable. 

SUMMARY OF THE INVENTION 

According to one embodiment of the present invention, an 
image interpolation method based on a sub-unity coordinate 
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shift technique is provided. In the approach, an original 
digital image in the spatial domain (e.g., Cartesian coordi- 
nate system (x,y)) is transformed into a frequency domain 
representation via a Fourier transform. The Fourier trans- 

s form representation is then modified by phase shift terms 
corresponding to image shifts in the spatial domain with 
sub -unity distances matching the locations where the image 
values need to be restored. The original and the shifted 
images are then interspersed together, yielding an interpo- 

10 lated image. The present approach can achieve an interpo- 
lation quality as good as that with the sine function since the 
sub-unity shift in Fourier domain is equivalent to shifting the 
sine function in spatial domain, while the efficiency, thanks 
to the fast Fourier transform (FFT), is very much improved. 

is In comparison to the conventional interpolation techniques 
such as linear or cubic B-spline interpolation, the interpo- 
lation accuracy is significantly enhanced. The method is 
especially applicable to two and three-dimensional com- 
puter tomography (CT) and magnetic resonance imaging 

20 (MRI) images, as well as other medical and non-medical 
digital images. 

It is therefore an advantage of the present invention to 
provide a method for interpolating digital images to recover 
lost image data due to digitization. 

25 It is a further advantage of this invention to provide a 
method for processing digital images that is computationally 
eflScient. 

It is yet a further advantage of the present invention to 
30 provide a method for processing digital images that is 
particularly suited to medical images such as, for example, 
computer tomography (CT) and magnetic resonance imag- 
ing (MRI). 

BRIEF DESCRIPTION OF THE DRAWINGS 

35 

In the accompanying drawings which are incorporated in 
and constitute a part of the specification, embodiments of the 
invention are illustrated, which, together with a general 
description of the invention given above, and the detailed 
40 description given below, serve to example the principles of 
this invention. 

FIG. 1 is a flowchart of the logic of the present invention 
that can be implemented on a computer system. 

FIG. 2 is a diagram illustrating the a 2 dimensional image 
45 undergoing sub -unity shift process of the present invention. 

FIG. 3 is a diagram illustrating the sub -unity process of 
the present invention in the context of a 3 dimensional 
image. 

50 DETAILED DESCRIPTION OF ILLUSTRATED 
EMBODIMENTS 

The detailed description shall be made by way of refer- 
ence to the drawing Figures and text herein. Prior to dis- 
ss cussing the drawing figures, the present discussion shall 
focus on the several important principles of the present 
invention. 

In this regard, a number of preliminary definitions shall be 
discussed. More specifically, the sets of integer numbers, 

60 real numbers, and complex numbers shall be denoted by Z, 
R and C respectively. Moreover, X <= Y shall be defined as 
set X is a subset of set Y and /:X-*Y shall be defined as a 
function f from a set X into a set Y. The Cartesian product, 
X 1 xX 2 x . . , xX„, of n independent sets shall also be defined 

65 as an n dimensional coordinate set. If all the sets are subsets 
of the set of real numbers, the product can be denoted in 
short by R". Given coordinate and value sets as subsets of 
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the set of real numbers, XcR" where ne{2,3}, and LcR, CPU running IRIX Version 63 in 128 MB RAM. As will be 

respectively, then an L valued image a on X is the graph of described, the present invention utilized several forward and 

a function a:X->L. Thus an L valued image a on X is of the backward Fourier transformations. The discrete implemen- 

form: . tation must be done with care to avoid artifacts due to signal 

5 discontinuities in both spatial and frequency domain. In 

0) some cases, it is helpful to window an image in spatial 

, x T T , . L / / • domain before the Fourier transformation if there are dis- 

where a(x)eL. In this graph, (x, a(x)) represents an image contimlities between ^ itc ^ of ^ ^ fleId 

element of the image a, where x is the pixel locauon and a(x) Referring now to FIG. 1, a flowchart 100 is shown 

the pucel value at location x. After d.gitizaUon, the union of w fllustrati thc lo ^ c of me , In particu]ari 

all the centers of.the samphng cells, referred to as a sampling ^ , ic cogences with step 102 where an original image 

grid, is usuaUy depicted by a set of integer numbers or rfd js ^ t m , 0 the tem ^ ■ fa rcferably m 

xeXe Z", while the _prxel values can be mteger, real, or ^ spatial domaiD) which can be deflned by a Cartesian 

complex, a(xj£Lc |Z,R,C). . coordinate or other conventional coordinate system. In step 

In this regard image interpolation is used to recover the „ 1M ^ final . resolution ^ entered whicb fa fer . 

missing image elements due to the process of digitization. If M ter ^ ^ ori inal ^ resohnion (e g > 2 ^ 

the discretization of the image field satisfies the Nyquest magnificaU o Dj etc .). step i 06> the interpolation ratio along 

samphng criterion the pixel value at any location inside the each coordinale de t ermine d. For example, to expand a 2-D 

imaging field can be accurately determined. Let a denote a { on a 512x sn grid (i.e., pixel or voxel) to that on a 

digital image acquired from an object or a scene. Given a set 2Q 4096x4096 rid) tbe interpola ti OD ratio k 8 in each direction . 

of parameters meM<=R , characterizing the coordinate To j,,,^,^ m6 missi sliccs of a set of CT ^ from 

translation, the new pixel coordinate is a function, ?5 slices tQ 3y5 slices> ^ interpolltioil ratio ^ 5 in the 

y:XxM-Y, where Y<=R . Each pixel location yeY can be longitudinal Omt&m and 1 within the slice. The interpola- 

calculated. t * on rat - Q determines how many shifted images need to be 

_ x m ^ 25 generated by the system. As for the 2-D example wilh an 

*~ X m interpolation ratio of 8x8, 63 shifted images are needed. 

The pixel value which has lost in the process of digitization In step 108, the original image is transformed, as 

can be approximated by the present invention. If image a is described above via a Fourier Transform, from the spatial 

finite and band limit, its Fourier transform ;f:XxS-»L, where domain into a frequency domain representation of the 

ScR M and LcC exists: 30 i ma g e * If tne image is not continuous from side to side, 

windowing of the image before the transformation may be 

F{sy\a{x)e- t2xa dx (3) required. In step 110, the Fourier transform of the original 

image is multiplied by a phase term t~ i2nms in the direction 

Where, F(s)eL and seS . The shift theorem of the Fourier i mage ^ to be shifted, where m is the reciprocal of the 

transformation states that if a(x) has a Fourier transform 35 interpolation ratio. 

F(s), then a(x-m) has the Fourier transform Q- £23ims F(s). On [ n step 112, an inverse Fourier transform is executed on 

the other hand, one can multiply the Fourier transform of the th e shifted frequency domain representation to generate a 

original image by a phase term e -1 " 2 ^. The inverse Fourier shifted spatial domain representation (i.e., image). The 

transform of the product yields the shifted image b:Y-*L « rea i" portion of the inverse transform is a shifted version of 

from the original image a: 4 q the original image. Steps 110 and 112 are repeated until all 

the images covering all the missing points are obtained, as 

H*W*~' F(s)e^ds (4) directed by step 114. For example, to interpolate- a 2-D 

Since o'^ F(s) is a real even and imaginary odd function, ^ from 25 f X ^ 6 S rids t0 three shifted 

the resulting image reduces to a real function b(x)eLc R. T^^F^T^ * T^l n n ?<l 

If a digital im/ge a, satisfies the Nyquest sampling crite- 45 <?' 5 > ( 00 > °- 5 ); and °5>- f the 

,? P . t c • + *u • i « have been generated, step 116 interlaces the elements ot all 

non and its Fourier transform exists, the image values at any " . 6 . ' v i , " , . 

locations other than the samphng points can be uniquely '^.images Rowing heir location indices^ The location 

determined with Equation (4). Suppose that N-l shifted Jf^fL"? de k fined h J. ** ** .direction and amount (i.e., 

images, a,, a 3 . . . .and a„, are obtained from a 1( retaining all ?f sh f t m *>c vertical direction etc.) The result is an 

the missmg image elements. Image interpolation involves to 50 jnlftpolited version of the ongmal image at the new reso- 

re-organizing the image elements from the N images. In this ^?JL e I e .: 1 . e , 4 . . , 4 

regard, the interpolated image can be denoted by h:X-L. , FI °- 2 d . hurt ^ ™. exam / le of h ™ ,0 ? te , M a 

& r 6 3 two-dimensional (2-D) image from an NxN grid to a 2Nx2N 

h~{(x,h(x)ys ex} (5) 8"^- As described above, the original image 202 is first 

55 transformed into a frequency domain representation in the 

where, the coordinate set is an ordered union of the coor- form of a Fourier transform. The Fourier transform is then 

dinate sets of. the N images, X={X 1 UX 2 UX 3 U. . . multiplied by a phase term e " t ^ t0 ' 5 ' r along the abscissa 

UX^JcR", with the image elements interspersed in direction only, the ordinate direction only, and both the 

sequence, and the value set becomes the union of the value abscissa and the ordinate directions respectively, as repre- 

sets, L={L 1 UL 2 UL 3 U. . . Ul^JcR. 60 sented by grids 204, 206, and 208. The inverse Fourier 

The above definitions describe a general approach to data transforms of these three functions are three shifted images, 

interpolation according to the present invention utilizing as represented by grids 204, 206, and 208, by 0.5 pixels in 

sub -unity coordinate shift methods. the horizontal and the vertical directions, and by 0.5^2 pixels 

The present invention was implemented with an IDL in the diagonal direction respectively. Re-arranging the 

graphics programming system from Research Systems Inc., 65 image elements yields an expanded or interpolated image, as 

of Boulder, Colo., on a SGI o2™ Workstation from Silicon shown by grid 210. In grid 210, the "o"s represent image 

Graphics Inc., of Mountain View, Calif., with a MIPS R5000 elements of the original image, "A"s denote elements shifted 
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by half a pixel of the original image in the horizontal 
direction, and " 0 "s denoted elements shifted by half a pixel 
in the vertical direction, and "+"s denote elements shifted by 
0.5^ pixels in the diagonal direction (i.e., 0.5 in the vertical 
and 0.5 in the horizontal direction). 

The present invention has been applied to X-ray computer 
tomography (CT) and magnetic resonance imaging (MRI), 
which can provide patient volumetric information in mul- 
tiple slices. It is frequently required to interpolate the image 
slices, so that the image elements are isotropic. As shown in 
FIG. 3, this can be done by shifting the image as a set of 3-D 
volume data in the longitudinal direction with sub -unity 
distances of, for example, 0.33 voxels, and then interlacing 
the slices of the original and all the shifted images. Since the 
present method is based on a sub-unity shift, such shifts can 
range from greater than 0 to anywhere less than one, 
depending on the number of coordinate shifts desired. 

While the present invention has been illustrated by the 
description of embodiments thereof, and while the embodi- 
ments have been described in considerable detail, it is not 
the intention of application to restrict or in any way limit the 
scope of the appended claims to such detail. Additional 
advantages and modifications will readily appear to those 
skilled in the art. For example, the present invention can be 
implemented on a general purpose Personal Computer or 
networked computer. Other possible medical images that 
can be analyzediiy the present invention include ultrasound 
images and histological images. Additionally, the present 
invention may take the form of a computer-readable medium 
that has the logic or instructions necessary for instructing a 
computer to execute commands and processes according the 
present invention. Such computer readable mediums 



10 



20 



25 



include, for example, magnetic disks, read-only memories 
(ROMs), programmable read-only memories (PROMs), 
electronically programmable read-only memories 
(EPROMs), electronically erasable programmable read-only 
memories, and compact-disk read-only memories (CD- 
ROMs). Therefore, the invention, in its broader aspects, is 
not limited to the specific details, the representative 
apparatus, and illustrative examples shown and described. 
Accordingly, departures may be made from such details 
without departing from the spirit or scope of the applicant's 
general inventive concept. 
We claim: 

1. A method of generating an interpolated image com- 
prising the steps of: 

(a) reading at least one original digital image in the form 
of an original spatial domain representation; 

(b) entering a final image resolution; 

(c) determining an interpolation ratio based on the origi- 
nal digital image and the final image resolution; 

(d) transforming the original spatial domain representa- 
tion into a frequency domain representation; 

(e) shifting the frequency domain representation by a shift 
value; 

(f) transforming the frequency domain representation to a 
shifted spatial domain representation; and 

(g) interlacing the original spatial domain representation 
and the shifted spatial domain representation to arrive 
at an interpolated image. 
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